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NUMERICAL DISPERSION ANALYSIS OF THE CONVECTED 
HELMHOLTZ EQUATION 


OHSUNG KWON AND IMBO SIM 


Abstract. We present the numerical dispersion effects in solving the con- 
vected Helmholtz equation by the conforming and nonconforming quadrilat¬ 
eral finite elements. Particularly, we evaluate the dispersion relations for the 
numerical schemes. The dispersive behaviors are analyzed by focusing on the 
Mach number and the angular frequency. Numerical experiments are con¬ 
ducted to verify the relations between the numerical dispersions and the com¬ 
putational errors. 


1. Introduction 

In this paper, we consider the convected Helmholtz equation in the presence of 
a mean flow. This equation is generated from the linearized Euler equations by 
reducing it for the pressure field to describe a propagating wave to the mean flow 
(see [16]). So it has been applied in sciences and engineering problems, for instance, 
in aeroacoustics (din])- In addition, various studies have been devoted as follows: 
Becache et al. used perfectly marched layers for the convected Helmholtz equation 
to design efficient numerical absorbing boundary conditions in |3]. Casenave et al. 
in [^ computed a linear acoustic wave propagation at a fixed frequency in the 
presence of flow using the coupled BEM-FEM. An algebraic subgrid scale finite 
element method was presented by Guasch et al. to improve the accuracy of the 
Galerkin finite element solution in [5]. However, issues concerning the dispersion 
properties of the finite element method remain unsolved yet. 

The dispersion relation concerns the angular frequency of a wave to the wavenum¬ 
ber and the Mach number. From this relation the phase and group velocities of the 
wave are derived. Sometimes they lead instabilities due to the opposite signs (see 
[nuttH]). Hence the dispersion analysis is an important issue. For the Helmholtz 
equation, there are numerous studies concerning the dispersion properties. Harari 
et al. in [siiinj used the Galerkin least squares to solve the Helmholtz equation and 
[Hill (HI developed the generalized finite element method. Other approaches are 
appending the element boundary residuals to the Galerkin approximation in |15j 
and the nonconforming element in |19| . 

The aim of this paper is to analyze the numerical dispersion behaviors of the 
convected Helmholtz equation using the conforming and nonconforming finite el¬ 
ement methods. Particularly, we investigate the difference of the continuous and 
numerical angular frequencies — and fine that it relates to the numerical errors. 
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2. Numerical dispersion relation of the convected Helmholtz 

EQUATION 

Let us consider the convected Helmholtz equation in a uniform mean flow on the 
domain H S 

—div(HVp) — 2ioj M • Vp — w^p = 0 in H, 

V ■ (AVp) — iujp = g on 

Here, M = (M, 0) for the Mach number M and H is a 2 x 2 diagonal matrix with 
Hii = 1 — and A22 = 1- Then the variational form of the problem (EH) is to 
find p S H^(H) such that 

(2.2) (HVp, Vu) — 2iu;(M • Vp, n) — a;^(p, n) — ia;(p, u) = ( 5 , i;), Vn £ H^(H), 

where (•,•) and (•,•) are the L^-inner products such that {u,v) = J^uvdxdy and 
(u,v) = fg^uvds, respectively. 

Next we will investigate the dispersion relations by the conforming and noncon¬ 
forming finite element methods using quadrilateral elements of the lowest order. 
In particular, two numerical schemes are used: PI conforming (Pl-C) method, 
Rannacher-Turek nonconforming (RT-NC) method. Here, there are two types of 
RT elements. We set RTl element by Rannacher-Turek element with the midpoints 
of edges as the degrees of freedom and RT2 element by one with mean integrals 
over edges (see [T7]b 



Figure 1. Computational region H = [—h,/i]^ 


We first evaluate the dispersion relation by RTl-NC method. To do so, let p^ 
be the numerical solution of the problem (12.21) such that 

(2.3) p^ = exp{i(kix + k2y)}, 

which is a plane wave propagating with the numerical wave vector )■ 

Since the plane waves have the same structure on each quadrilateral mesh, we 
restrict the domain fi = [—h, as Figure [U Then H contains 12 midpoints ruj 
and 4 rectangles Bi. The central point O is not used in RTl-NC method, but it 
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needs to define a global test function ipc- In particular, ipc is defined by 

Vb\Ri + on TOi, 

Vt\r 2 + Vi\rz on TO 2 , 

<Pb\Ri + on TO 3 , 

‘f>r\R^ + Vl\Ri on TO 4 , 

where (pj{j = r,l,t,b) is a local basis function on Ri and the subscript r, I, t, and 
b mean the right, left, top, and bottom on i?/, respectively. Similarly, we denote 
(fij the basis function at rrij for j = for examples, ipi = if>h\Ri + 

Using these basis functions Lpj of RTl-NC method, the numerical solution is 
represented by ^or constants pj. Inserting p^ and po into (12.21) 

and using the middle point rule to approximate the inner products on the boundary, 
we have 

12 

^Pj{{AVpj,ypG) - 2za;(M • Vpj,pG) - = 0. 

i=i 

By the direct calculation of the inner products, we have 

,,2l2 4 12 

j=l j=5 

- -pi) + 2(pio -Pe) + 2 (p 9 -ps) + (ps -P12) + (pii - Pv)} 

4 12 

+ ~'^Pj + M'^{p 5 +P 6 +P 9 +P10 - 2 (p 2 +P4)} = 0 . 

i=i j=5 

By (ESI), PjPj = exp{i{kix + kl^y)}. Using the properties pj{mj) = 1 and 
Pj{mi) = 0 for j ^ I, Pj is expressed in terms of kj and h, e.g., pi = exp{—ikih/2). 
So, letting Cj = cos{kjh/2) and Sj = sm{kjh/2) for j = 1, 2, the dispersion relation 
of RTl-NC method is given by 



(2.4) 

where 


4M/ 


/ 3(C'i+C2)(1-CiC2)-M252C'2 

1 2 M2(Ci+C2)(2 + CiC2) 


G, = 


^iG 2 ( 2 Gi+G 2 ) 

'' " (Gi+G2)(GiG2+2)- 

Let Cj = cos{kjh), Sj = sm{k^h) for j = 1,2. By the same method used in (12.4L 
we obtain the dispersion relations for RT2-NC and Pl-C methods such that 
(2.5) 


OJ = 


6M 

ST 

3M 


w =- 


G 2 + 


G 3 + 




^ 2 fc(‘52(l-Gi) + fc^^i(l-G2)-M2fc(^52(l-Gi) 
" 3 M^{k'lS2{5 + Ci) + k^Si{5 + C2)) 


\ 


2 2 2 G 1 G 2 + Gi + G 2 - 4 + M2(2 - G 1 G 2 - 2Gi + Ga) 


G|- , 
3 3 


M 2 (CiG 2 + 2 (Gi+G 2 ) + 4) 
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M=0.3, 6=0 M=0.3, 6=0.25 7 : 





M=0.9, 6=0 


M=0.9, 6=0.25 7: 




Figure 2. Dispersion quotient bias H for M = 0.3,0.6,0.9 and 
9 = 0,0.2571 


where 

G fci §-^2 + ^(1-C0(l + C2) ^ gi(2 + C2[ 

k'lS2{5 + Cl) + kl^Siib + C 2 )' ^ CiC2 + 2(Ci+C2)+4’ 

respectively. 

Remark 2.1. The dispersion relation of the continuous problem (EB is of the form: 

(2.6) oj = kiM + |k|, 

which is derived from the characteristic equation of EB by ESD- 

3. Analysis of the numerical dispersion relation 

To derive the dispersive behaviors of (I2.4I1 - (I2.5I1 . we define a new variable 9 that 
is the angle for the direction of the wave propagation so that k = A:(cosd, sind) 
with k = |k|. Then (12.611 is represented by w = k{l + M cos 9). As a result, the 
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M=0.3, e =0.7 5k 



M=0.3, 6 = k 



— RTl-NC 

- - RT2-NC 

■ ■ Pl-C 





0.0 H 0.3 

M=0.6, e = K 



M=0.9, 9 = k 



Figure 3. Dispersion quotient qp bias H for M = 0.3,0.6,0.9 and 
9 = O.TStt, tt 


dispersion relations in (I2.4I) - (I2.5I) are reformulated in terms of Then we can 
expand w in terms of h such that 

OO 

(3.1) uj = + ^ Aj{M, 9) 

i=i 

where Aj{M, 9) is a function depending on M and 9. This yields that converges 
to the continuous angular frequency w as h —> 0. 

To compare the dispersive behaviors of the numerical schemes, we dehne two 
dispersion quotients in [^: 

qp = Vp/v^, qg = Vg/v^, 

where Vp and Vg are the phase and group velocities in the direction of the numerical 
wave vector such that Vp = w/kA and Vg = |9a;/9k^|, respectively. Similarly, 
the numerical phase and group velocities are = uAjkA and 
respectively. Then qp and qg measure the errors in the phase and group velocities. 




















6 


OHSUNG KWON AND IMBO SIM 


M=0.3, 9=0 


M=0.3, 9=0.25K 




M=0.6, 9=0 


M=0.6, 9=0.257: 




M=0.9, 9=0 


M=0.9, 9 =0.25K 




Figure 4. Dispersion quotient bias H for M = 0.3,0.6,0.9 and 
9 = 0,0.2571 


Using (13.IL we have qp = 1 + . To observe the effects of 

M and 9 clearly, we set H := uj^h G [0,0.3], since the numerical experiments will 
be done on this interval. Since the mean flow passes horizontally, the dispersion 
quotient qp is symmetric with respect to x axis, so it suffices to consider 9 G [0, tt]. In 
particular, we use four angles of 9 such as 0, 0.257r, 0 . 75 tt, and tt, since the different 
behaviors of the numerical schemes are well observed on these angles. For Mach 
number M, we set M = 0.3,0.6, and 0.9. The asymptotic behavior of qg is also 
similar to so we use the same setting. Under these conditions, the dispersion 
quotients qp and qg are illustrated as Figures [2HS1 respectively. 

Figures [2]in] show the effects of M and 9 to the dispersion quotients of numerical 
schemes. For RT-NC method, the dispersion quotients get away to 1 as M and 
9 are larger. Specially, it is significantly large at M = 0.9 and 9 = tt. For Pl-C 
method, it shows the anisotropy behavior to M. For large 0, it appears the similar 
behaviors to RT-NC method. Meanwhile, the dispersion quotients approach to 1 
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M=0.3, 6=0.7571 






0.0 H 0.3 0.0 

M=0.6, 6 =0.7 5k 


H 0.3 

M=0.6, e = K 




Figure 5. Dispersion quotient qg bias H for M = 0.3,0.6,0.9 and 
9 = O.TStt, tt 


as M is larger for small 0. This means that Pl-C method is more efficient than 
RT-NC method for small 9 and large M, though it does not for other cases. 

In fact, the behaviors of the dispersion quotients are determined by the differ¬ 
ence of the numerical and continuous angular frequencies. In next theorem, we 
investigate the behavior of the difference — uj using techniques in mils]. 

Theorem 3.1. Let oj and uj^ be continuous and numerical angular frequencies, re¬ 
spectively. We define the dispersion error by \(jj^—uj\. Then \uj^—uj\ = Ai{M, 9)uj^h^-\- 
0{uj^h^), where 

|2(1 -I- cos(40)) -I- 4M(cos(30) — 3cosd) -I- M^(cos(40) — 6cos(20) — 7)| 

384(1-hM cos 6»)3 ’ 

|2(1 -I- cos(40)) -I- 4M(cos(30) — cos0) -I- M^(cos(46*) — 2cos(20) — 3)| 

^ “ 192(1-hM cos 0)3 ’ 

|3 -I- cos(40) — 4M^ cos^ 9\ 


96(1 -|- M cos 0)3 
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6=0 9 =0.25 K 



6=0.757: 




0.3 0.6 0.9 

M 


Figure 6. Graphes of Ai{M,0) bias M ioi 9 = 0,0.257r,O.TS tt,tt 


which are the leading coejficient functions for RTl-NC, RT2-NC, and Pl-C meth¬ 
ods, respectively. 

Generally, if h is sufficiently small so that uJ^h? < 1 for large w, then the main 
behavior of the dispersion error is determined by Ai {M, 6). This fact is confirmed 
by that the effects of M and 6 to the dispersion quotients in Figures [5]-[S] follows 
the behaviors of Ai{M, 9) in Figure [51 

To check the effects of the dispersion errors to the computational errors of the 
numerical schemes, we define the norm: 

(3-2) IIIpIII := II V^Vp|||2(f2) + w^||p|||2(f2), 

which is the energy norm of (12.21) for the real part. To calculate the numerical 
solution, we need the boundary date g. It could be verified by = exp{i{kiX + 
kfy)} such that 

g = i{kv ■ (Aer) — ui) exp(zfcx • e^), 
where e^. = (cos6>,sin0) and x = {x,y). 

In order to observe the effect of Ai{M, 9) clearly, we set = I and solve the 
problem by growing w from 10 to 80. Figures |3 and [8] illustrate the numerical error 
Err := |||p — P^||| for M and 9, which yields that the error behaviors nearly follow 
Ai{M,9) in Figure [5] 

4. Comparison to the Helmholtz formulation 

The convected Helmholtz equation in (12.11) could be reformulated as the Helmholtz 
equation by the following lemma. 
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M=0.3, 6=0.25K 
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80 


M=0.9, 6=0 






Figure 7. Numerical errors bias w for M = 0.3,0.6,0.9 and 0 = 0, 0.257r 


Lemma 4.1. Let a{x) = . If we set u{x,y) := p{x,y)a{x), then 

u{x,y) is a solution of the problem: 

(4.1) — div(ylVM) — ® in LI. 

Proof. Let d = 1 — M^. Inserting u{x,y) = p(x, y)a{x) into (14.11) . we have 
^2 

-div(AVu) - = -d{pa)xx - {poi)yy - uj'^pa/d 

= -d{pxxa + 2pxax + paxx) - PyyOi - uPpajd 

= —da(pxx + 2iwMpx/d — M'^uj^p/d^) — PyyU — ui^pa/d 
= a[- dpxx - Pyy - 2iuiMpx - UJ^p). 

□ 

Using Lemma 14.11 we can find the solution p of the problem (12.11) without con¬ 
sidering the convection term by solving gH). However the problem dni has a 
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M=0.3, 6=0.75?: 

400 

Err 

200 

. :., . 

!0 0 ) 80 


M=0.3, 9 = 71 



M=0.6, 6=0.7571 



10 CO 

M=0.9, 9=0.7571 
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80 


Figure 8. Numerical errors bias w for M = 0.3,0.6,0.9 and 0 = 0.757r,7r 


stability deterioration for large Mach number by the term 1/(1 — M^) in m- 
This phenomena is generated by the dispersion error — w as follows: 

We first calculate the dispersion relations of the problem dSD. (j2.3p 3jnd 

( 12 : 61 ) . the numerical solution of u is of form = exp{z(A: 3 a; + kl^y)} with 
fcg = (fc/ + fc^M)/(l — M^). Using the process stated in Section 2, the dispersion 
relations for RTl-NC, RT2-NC, and Pl-C methods are given by 


LO = < 



on AfU ~ - 1) 

^ ^ (C2 + Cs) + (2 + C2C3) 

^2^3(1 - C 2 )+ k3S2{l - Cs)^- M^k3S2il - C 3 ) 

k2S3{C2 + 5) + k3S2{C3 + 5) 

^g.^ _ ,^^ 4 - 2C2C3 -C2-C3+ MHC2C3 -C2 + 2C3 - 2 ) 

C 2 C 3 + 2 (C 2 + C 3 ) + 4 
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where C 3 = cos{k^h/2) and C 3 = cos{k^h). Under the same assumption to h and uj 
in Section 3, the main behavior of the dispersion error is determined by the leading 
coefficient Ai{M,9), so it suffices to consider Ai{M,0) for each numerical scheme. 
For RTl-NC method, Ai{M,9) is given by 


Ai{M,9) 


(4 + lOM^ + 28M4 - 7M^) + 4M(4 + IIM^ - M*) cos9 
768(1-M2)2(1 + Mcos6»)4 

4M2(11 - 6M2 + 2M4) cos(26I) + 4M(4 - 3M^ + M^) cos(36») 
768(1-M2)2(1 + Mcos6»)4 
(4 - 6M2 + 4M4 - M^) cos(46») 

768(1-M2)2(1 + m cos 6»)4 ' 


For RT2-NC method, it is the double of one for RTl-NC method. For Pl-C method, 
it is 

I cos'* 9 + AM cos^ 9 + 6 M^ cos^ 9 + AM^ cos 9 + -|- (1 — sin^ 9\ 

^ ~ 24(1-M2)2(1-pM cos 0)4 ■ 


All Ai{M, 9) blow up as M —>■ 1, so it leads to instabilities for large M. 


5. Conclusion 

We have analyzed the dispersive behaviors of the convected Helmholtz equation 
by the conforming and nonconforming finite element methods. Particularly, the 
dispersion relations of the numerical schemes are derived on the quadrilateral mesh 
and it is shown that the numerical dispersions converge to the continuous ones as 
/i —>■ 0. We also observed the effects of the dispersion error — w| to the dispersion 
quotients and the numerical errors. 

Consequently, our result informs the effects of the Mach number M and the 
angular frequency oj to the numerical errors. So it provides the guidelines for the 
selection of an appropriate mesh size in terms of M and w in solving numerically. 
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